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ABSTRACT 

Over the past two decades, helioseismology has revolutionized our understanding of the interior structure and 
dynamics of the Sun. Asteroseismology will soon place this knowledge into a broader context by providing 
structural data for hundreds of Sun-like stars. Solar-like oscillations have already been detected from the ground 
in several stars, and NASA's Kepler mission is poised to unleash a flood of stellar pulsation data. Deriving 
reliable asteroseismic information from these observations demands a significant improvement in our analysis 
methods. In this paper we report the initial results of our efforts to develop an objective stellar model-fitting 
pipeline for asteroseismic data. The cornerstone of our automated approach is an optimization method using 
a parallel genetic algorithm. We describe the details of the pipeline and we present the initial application to 
Sun-as-a-star data, yielding an optimal model that accurately reproduces the known solar properties. 

Subject headings: methods: numerical — stars: interiors — stars: oscillations 



1. MOTIVATION 

Most of what we can learn about stars comes from obser- 
vations of their outermost surface layers. We are left to infer 
the properties of the interior based on our best current un- 
derstanding of the constitutive physics. The exception to this 
general rule arises from observations of pulsating stars, where 
seismic waves probe deep through the interior and bring in- 
formation to the surface in the form of light and radial veloc- 
ity variations. The most dramatic example is the Sun, where 
such observations have led to the identification of millions of 
unique pulsation modes, each sampling the solar interior in a 
slightly different and complementary way. The radial profile 
of the sound speed inferred from inversion of these data led 
to such precise constraints on the standard solar model that, 
before the recent controversy over heav y element abundances 
(lAsplund et al.l2004l:ICaffau et alJ2008l) . the observations and 
theory agreed to better than a few parts per thousand over 90% 
of the solar radius (Christensen-Dalsgaard 2002). Although 
such detailed inversions are not currently possible for other 
stars, pulsation data do allow us to determine the global prop- 
erties and to probe the gross internal composition and struc- 
ture, providing valuable independent tests of stellar evolution 
theory. The reason for this qualitative difference is simple: 
if we could move the Sun to the distance of the nearest star, 
many of the pulsation modes would no longer be detectable. 
Without spatial resolution, only those modes with the lowest 
spherical degree (Z < 3) would produce significant variations 
in the total integrated light or the spectral line profiles^. These 
are also the modes that probe deepest into the stellar interior, 
collectively sampling the physical conditions from the core to 
the photosphere. 

The solar oscillations are excited by the acoustic noise 
generated by the near-surface convection, with nearly sonic 
speed. Thus we expect to find similar oscillations in all stars 
with vigorous outer convection zones. A first unambiguous 
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detectio n of solar-like o scillations in another star was in fact 
made by lKjeldsen et al.l ([1995) . Recent improvements in our 
ability to make high-precision radial velocity measurements 
from the ground have been driven largely by efforts to detect 
extra-solar planets. This has led to the detection of solar-like 
oscillations in more than a dozen main-sequence and subgiant 
stars; in addition, large-scale ground-based photometric cam- 
paigns have shown evidence for solar-li ke oscillations in nu- 
mero us giants (for a recent review see iBedding & KjeldsenI 
l2007h . The oscillation amplitudes and the frequency of maxi- 
mum power in these stars agree reasonably well with our the- 
oretical expectations — ^but in some cases the mo de hfetimes 
appea r to be significantly shorter than expected (IStello et al.l 
l2004l) . suggesting that our knowledge of the convective driv- 
ing and damping physics is incomplete. Scintillation in the 
Earth's atmosphere severely limits our ability to detect the 
parts-per-million light variations due to these pulsations in 
main-sequence stars. However, space-based photometric pro- 
grams designed to detect extra-solar planet transits also have 
the sensitivity to document the s tellar pulsatio n signals. Thus, 
the French-led CoRoT mission (iBaglin et al.l l2006) is provid- 
ing extensive data of very high quality on oscillations in a 
broad range of stars. The field is progressing very rapidly and 
many new observations will become available in the next few 
years, particularly after the launch of NASA's Kepler mission 
in 2009. With this in mind, it will be a distinct advantage to 
have in place the computational methods that will allow us to 
maximize the science return of these data. This will lead us 
quickly to a deeper understanding of the solar oscillations in 
the context of similar pulsations in other stars, and will yield 
unprecedented constraints on the formation and evolution of 
stellar and planetary systems. 

The likely excitation mechanism by turbulent convection 
near the surface creates a broad envelope of power with 
a peak that scales with the acoustic cutoff frequency (see 
iBrown & GillilandllI994l) . Within this envelope a large frac- 
tion of the predicted low-degree pulsation modes are excited 
to detectable amplitudes, leading to readily identifiable pat- 
terns of peaks in the power spectrum characterized by the 
large separation {/S.vi = Vnj — u„-\ i), and the small separa- 
tion {5vi^i+2 = 1^,1.1-1^,1-1,1+2), where n is the radial order and 
/ is the spherical degree of the modes with oscillation fre- 
quency Vnj. Without any detailed modeling, these overall 
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patterns immediately lead to an estimate of the mean density 
of the star and can indicate the presence of interior chemi- 
cal gradients that are a good proxy for the stellar age. But 
a full analysis must include a detailed comparison of the in- 
dividual frequencies with theoretical models. One compli- 
cation with such a comparison is the existence of so-called 
surface effects, which appear as systematic differences be- 
tween the observed and calculated oscillation frequencies that 
grow larger towards the acoustic cutoff frequency. Surface ef- 
fects arise primarily due to incomplete modeling of the near- 
surface layers of the star where convection p lays a major role 
(IChristensen-Dalsgaard & Thompsonlll997h . Addressing this 
inherent deficiency in our ID models would require (among 
other things) that we substitute the results of extensive 3D 
calculations for the parameterized mixing-length treatment of 
convection that is currently used in nearly all stellar evolu- 
tion codes. Alternately, we can make an empirical correctio n 
to the calculated frequencies following^ Kieldsen et al.l (l2008l) . 
who recently devised a method for calibrating the functional 
form of surface effects using solar data, and then obtaining 
corrections for other stars based on the frequencies of a ref- 
erence model, scaled to correct for the different mean density 
(see AppendixlAl. 

We have developed an objective and automated method of 
fitting stellar models to the asteroseismic data soon expected 
to emerge from NASA's Kepler mission. This will lead to re- 
liable determinations of stellar radii to help characterize the 
extra-solar planetary systems discovered by the mission, and 
stellar ages to reveal how such systems evolve over time. For 
the asteroseismic targets that do not contain planetary com- 
panions it will allow a uniform determination of fundamental 
physical properties for hundreds of solar-type stars, providing 
a new window on stellar structure and evolution. In ^we 
provide an overview of the stellar models as well as the opti- 
mization and analysis techniques that establish the foundation 
of our fitting method. We describe our initial numerical exper- 
iments to calibrate the method with Sun-as-a-star data in ^ 
and we outline our plans for further validation in ^ 

2. COMPUTATIONAL METHOD 

The Kepler mission will soon yield precise high-cadence 
time-series photometry of hund reds of pulsating stars every 
few m onths for at least 3.5 years (IChristensen-Dalsgaard et all 
l2007h . We will then face the challenge of determining the fun- 
damental properties of these stars from the data, by attempt- 
ing to match them with the output of computer models. The 
traditional approach to this task is to make informed guesses 
for each of the model parameters, and then to adjust them it- 
eratively until an adequate match is found. The volume of 
asteroseismic data that will emerge from the Kepler mission 
calls for a more automated approach to modeling that initially 
explores a broad range of model parameters in an objective 
manner The cornerstone of our model-fitting approach is 
a global optimization method using a parallel genetic algo- 
rithm. The result of the global search provides the starting 
point for a local analysis using a Levenberg-Marquardt algo- 
rithm with Singular Value Decomposition, which also allows 
us to explore the information content of the observables and 
the impact of including additional obs ervational constraints 
(iBrown et alJII994tlCreevey et al.ll2007h . 

2.1. Stellar Models 

We hav e recently adapted the Aarhus stellar evolution code 
(ASTEC: IChristensen-Dalsgaardll2008ah and adiabatic pulsa- 



tion code ( ADIPLS : fChristensen-Dals gaardi l2008bh to inter- 
face with the parallel genetic algorithm. These are essen- 
tially the same models that were developed for the analy- 
sis of helioseismic data, and are t he source of Model S of 
Christ ensen-Dalsgaard et aP (Il996h . which has been used ex- 
tensively as a reference model for solar inversions. Using 
these models for the analysis of pulsations in solar-type stars 
will provide some internal consistency in our understanding 
of solar-like oscillations. Briefly, t hese stellar mode ls use 
the OPAL 2005 equation of state (see Rogers et al. 1996) and 
the most recent OPAL opacities (see Iglesias & Rogers 1996), 
supplemented by Kurucz opacities a t low temperatures. The 
nuclear reaction rates come from iBahcall & Pinsonneaulll 
( 1995) , convection is de scribed by the mixing-length theory 
of Bo hm-Vitensd (|1958), and we included the effects of he- 
lium settling as described by Michaud & Proffit ( 1993). 

Each model evaluation involves the computation of a stel- 
lar evolution track from the zero-age main sequence through 
a mass-dependent number of internal time steps, terminating 
prior to the beginning of the red-giant stage. Rather than 
calculate the pulsation frequencies for each of the 200-300 
models along the track, we exploit the fact that the aver- 
age frequency spacing of consecutive radial overtones {/S.vq) 
in most cases is a mono tonically decreasing function of age 
jChristensen-Dalsgaardl 1 1 9931) . Once the evolution track is 
complete, we start with a pulsation analysis of the model at 
the middle time step and then use a binary decision tree — 
comparing the observed and calculated values of (Ai'o) — to 
select older or younger models along the track. In practice, 
this recipe allows us to interpolate the age between the two 
nearest time steps by running the pulsation code on just 8 
models from each stellar evolution track. 

2.2. Global Search 

Since we are interested in developing a general-purpose 
modeling tool for asteroseismic data from the Kepler mission, 
we need to select a global method for optimizing the match 
between our model output and the available observations of 
any given star Using only observations and the constitutive 
physics of the model to restrict the range of p ossible values 
for each parameter, a genetic algorit hm (GA; ICharbonneaul 
ll995l:lMetcalfe & Charl3onneaul2003h can provide a relatively 
efficient means of searching globally for the optimal model. 
Although it is more difficult for a GA to find precise values 
for the optimal set of parameters efficiently, it is well suited 
to search for the region of parameter space that contains the 
global minimum. In this sense, the GA is an objective means 
of obtaining a good first guess for a more traditional local 
analysis method, which can narrow in on the precise values 
a nd uncertainties of the optimal mo del parameters. 

iMetcalfe & Charbonneaul (l2003h developed a fully parallel 
and distributed implementation of the PIKAIA genet i c algo - 
rithm^ that was originally written by ' Charbonneaul (Il995h . 
Metcalfe et al. (2000) used this modeling tool in the context 
of white dwarf asteroseismology, which ultimately led to a 
number of interesting physical results, including: (1) a pre- 
cise estimate of the astrop hysically impo rtant ('^C + "^He — + 
'^O) nuclear reaction rate (Metcalfe!l 2003l) . (2) the first unam- 
biguous detection of a crystallize d co re in a massive pulsat- 
ing white dwarf (iMetcalfe et alJl2004l) . and (3) asteroseismic 
confirmation of a key prediction of diffusion theory in white 

^ The parallel ve rsion of the PlKAlA genetic algorithm is available at 
|http://www.cisl.ucar.e dii/css/staff/travis/mpik:aia/ 
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dwarf envelopes dMetcalf j|2005h . The impact of this method 
on the analysis of pulsating white dwarfs suggests that seis- 
mological modeling of other types of stars could also benefit 
from this approach. 

Our implementation of the GA optimizes four adjustable 
model parameters; these are the stellar mass (M<,) from 0.75 
to 1.75 Mq, the metallicity (Z) from 0.002 to 0.05 (equally 
spaced in logZ), the initial helium mass fraction (Yq) from 
0.22 to 0.32, and the mixing-length parameter (a) from 1 to 
3. The stellar age (r) is optimized internally during each 
model evaluation by matching the observed value of (Ai'o) 
(see ^2. lb . The GA uses two-digit decimal encoding, so there 
are 100 possible values for each parameter within the ranges 
specified above. Each run of the GA evolves a population of 
128 models through 200 generations to find the optimal set of 
parameters, and we execute 4 independent runs with different 
random initialization to ensure that the best model identified 
is truly the global solution. This method requires about 10^ 
model evaluations, compared to 10'^ models for a complete 
grid at the same sampling density, making the GA nearly 1000 
times more efficient than a complete grid (currently 1 week 
of computing time, compared to many years for a grid). Of 
course, a grid could in principle be applied to hundreds of ob- 
servational data sets without calculating additional models — 
but the GA approach also gives us the flexibility to improve 
the physical ingredients in the future, while the physics of a 
grid would be fixed. 

2.3. Local Analysis 

Once the GA brings us close enough to the global solution, 
we can switch to a local optimization method. We imple- 
ment a modified Levenberg-Marquardt (LM) algorithm that 
uses Singular Value Decomposition (SVD) on the calculated 
design matrices to filter the least important information from 
the observables (some of which may be dominated by noise). 
This then provides an effective local inversion technique. LM 
is relatively fast and stable, and convergence typically occurs 
within 3-4 iterations. 

We treat the local analysis as a x^-minimization problem, 
where one seeks to find the set of parameters P that minimizes 




Here O,- and cr, are the / = 1 , 2, . . . ,M measurements and errors, 
while the C,- are the calculated model observables resulting 
from P. The LMh-S VD method requires an initial guess of P, 
and these are taken to be the results from the global search by 
the GA. The LMh-SVD analysis subsequently uses derivative 
information from the model at the current parameter values to 
calculate suggested parameter changes SP that will bring the 
model observables closer to the observations. 

We have three main motivations for implementing a local 
optimization method at the end of the global search. First, the 
GA has a limited resolution for each parameter, and the values 
that match the observations best are most likely between the 
fixed sample points. For example, the parameter Z near solar 
values is tuned by the GA along sample points spaced 0.0006 
apart, corresponding to a precision of roughly 3%. Consider- 
ing that the search resolution is limited for all of the model 
parameters, this significantly limits the precision of the global 
optimization. The resolution of the local analysis is limited 
only by the precision of the stellar evolution and pulsation 
codes, so we use it to adjust the models below the resolution 
of the GA search. 



Our second motivation for the local analysis is to quan- 
tify the final parameter uncertainties and correlations, and to 
probe the information content of the observables. We do this 
by calculating the derivatives of the model observables with 
respect to each of the fitted parameters, and then dividing each 
vector by the corresponding measurement error. This matrix 
can be referred to simply as the design matrix D. We subse- 
quently calculate the singular value decomposition of D for 
some of its very useful properties (see Appendix iBt. 

The final motivation for implementing a local analysis is to 
explore the effects of usin g different physical descriptions of 
the stellar interior (ICreeve y 2008, 2009). When the changes 
to the underlying physics are relatively subtle, we can assume 
that the global search by the GA using one set of assumptions 
will also provide a good starting point for a local analysis un- 
der slightly perturbed condit ions. By applying the techniques 
explained in ICreevevI (l2009l) . such an analysis could demon- 
strate that the observational data contain enough information 
to distinguish clearly between different choices for the equa- 
tion of state, for example. Ultimately, this technique could 
reveal discrepancies in the observables that indicate which 
physical description is most suitable for the star under inves- 
tigation. 

We apply the local analysis to each of the sets of optimal 
parameters identified by the four independent runs of the GA. 
The local analysis involves the gradual refinement of the op- 
timal parameters through the iterative application of several 
steps: 

1. Scan the Z parameter within ±0.002 of its original 
value, and perform a new minimization for each start- 
ing value of Z. We retain the set of parameters Pnew that 
result in the lowest value of the reduced x^- 

2. Starting with Pnew, scan the parameter within ±0.01 
of its initial value, and perform a new minimization for 
each starting value of Yq. Again we retain the Pnew that 
yields the lowest reduced x^- 

3. Rescan the Z parameter beginning with Pnew- We adopt 
as the optimal model the set of parameters that yield the 
lowest reduced from this final iteration. 

The reason for scanning Fq and Z becomes clear when we be- 
gin to understand the intrinsic parameter correlations, which 
are enhanced by the limited set of observations (see ^3.2\ . 

In most cases the best model found by the global search 
leads to the final best model after the local optimization. How- 
ever, sometimes a model that appears to be marginally worse 
at the end of the global search is improved more substantially 
during the local analysis. We take our final solution to be the 
best match from the four independent analyses, and for clarity 
we report the results of the corresponding global search even 
when it is not the best of the four models identified by the GA 
(cf. Table [T]i. The final uncertainties in the parameter values 
(M^,,Z,Yo,a,T) and in the model observables (reff,L*,7?,t) are 
calculated using SVD (see Appendix iBt. 

3. MODEL-FITTING EXPERIMENTS 

The overall goal of our model-fitting pipeline is to take a 
range of oscillation frequencies and other constraints as input, 
to identify and refine the model that best matches these obser- 
vations, and to produce the optimal values of several parame- 
ters and other characteristics of the model as output. To ensure 
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that our pipeline yields reliable results, we must begin by ap- 
plying it to data where the basic stellar properties are already 
known. The most obvious choice is the Sun, where high- 
quality Sun-as-a-star observations are available from multiple 
experiments. If we feed the pipeline a stellar-like set of solar 
data, we can judge the experiment a success if the pipeline re- 
turns the known solar properties within acceptable tolerances. 
To reach this goal, we must first optimize the efficiency of the 
search method by passing synthetic data through the model- 
fitting procedure 03.1l l, calibrate the differences between our 
models and real observations due to near-surface effects and 
ensure that the resulting empirical correction does not intro- 
duce any large systematic errors in the final model parame- 
ters 03. 2b . and finally quantify any differences in the derived 
model parameters due to the source and error properties of the 
solar data 03.3I ). The results of these experiments are listed in 
Table [T] and described in the following subsections. We will 
eventually want to validate this model-fitting pipeline using 
other stars that differ from the Sun, where the physical proper- 
ties are known with lower precision. However, because of the 
computation-intensive nature of our method (which currently 
requires about 1 week on 512 processors), we will consider 
here only the initial validation using solar data. 

3.1. Efficiency of the Search 

Although genetic algorithms are often more efficient than 
other comparably global optimization methods, they are still 
quite demanding computationally. Fortunately, the procedure 
is inherently parallelizable; we need to calculate many mod- 
els, and each one of them is independent of the others. So 
the number of available processors determines the number of 
models that can be calculated in parallel. Also, there is very 
little communication overhead; parameter values are sent to 
each processor, and they return either a list of observables or 
just a goodness-of-fit measure if the predictions have already 
been compared to the observations. The parallel version of 
the PIKAIA genetic algorithm is perfectly general, and did 
not require any structural modifications to interface with our 
stellar evolution and pulsation codes. 

The efficiency of genetic-algorithm-based optimization can 
be defined as the number of model evaluations required to 
yield the global solution, relative to the number of models that 
would be required for a complete grid at the same sampling 
density. In practice, a GA is usually hundreds or even thou- 
sands of times more efficient than a complete grid, and its per- 
formance is fairly insensitive to the few internal parameters 
that control its operation. We initially set these internal pa- 
rameters (e.g. population size, run length, crossover and mu- 
tation rates) based on our experience with white dwarf mod- 
els, but we also ran synthetic data through the optimization 
procedure in a series of "Hare & Hound" (H&H) exercises to 
ensure that the input parameters were recovered faithfully. 

The basic procedure for an H&H exercise is for one team 
member to calculate the theoretical oscillation frequencies 
and other observables for a specific set of model parameters 
(M,^,Z,lo,a,T). A subset of the predictions, typical of what- 
ever is available from actual observations, is then given to an- 
other team member who does not know the source parame- 
ters. The data are passed through the complete optimization 
method in an attempt to recover the source parameters without 
any additional information. The results of such exercises are 
used to quantify the success rate of the optimization method 
(the fraction of independent runs that lead to the known source 
parameters), and to improve its efficiency (minimize the num- 
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Fig. 1 . — The reduced of the best solution in the population as a function 
of generation number for four independent runs of the genetic algorithm on 
the H&Hl data set. All four runs converge to the global solution within 150 
generations. 



ber of model evaluations) if possible. 

For our first H&H exercise (with source parameters listed 
in the row labeled "H&Hl" in Table [T]i, we assume that typ- 
ical asteroseismic data from the Kepler mission will include 
twelve frequencies for each of the radial (/ = 0), dipole (/ = 1), 
and quadrupole (/ = 2) modes, with consecutive radial orders 
in the range n = 14-25^. Thus, we allowed the GA to fit a 
total of 36 oscillation frequencies. We assigned statistical un- 
certainties to each frequency by scali ng up the errors on the 
corresponding modes in BiSON data jChaplin et al.ll 19991) by 
a factor of 10, which is roughly what we expect from Ke- 
pler data (ffi, ^ 0.\ /iHz). We complemented this synthetic 
asteroseismic information with artificial data on the effective 
temperature and luminosity, with errors comparable to what is 
expec ted for stars in the Kepler Input Catalog ( Latha m et al.l 
l2Q0l reff = 5777±100 K, L*/L0 = l.OOiO.l)'^. This is the 
simplest conceivable test of the pipeline — if we generate a set 
of model data and then use the optimization method to match 
those data using exactly the same models, how quickly will 
the GA find the true values of the model parameters? The 
results of such a test verify the basic functionality of the al- 
gorithm, and tell us the approximate number of iterations (or 
"generations" of the GA) that we need to execute before stop- 
ping the search. 

The results of H&Hl for four independent runs of the GA 
are illustrated in Figure[T] Each convergence curve shows the 
reduced (hereafter x| = /{M-N - 1) where M is the 
number of constraints and is the number of parameters) for 
the best model in the GA population as a function of genera- 
tion number. All of the runs converge in less than ^-^150 gen- 
erations, in each case leading to the same mass as the source 

Note that only the identification of / is necessary for the operation of our 
pipeline. Information about the value of n is only provided for completeness. 

' Note that the stellar luminosity is a derived quantity based on the ob- 
served parallax and apparent magnitude. 
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model and values for the other parameters offset by less than 
~ 1 % from the source values (the best model from the GA is 
listed in the row labeled "global" under H&Hl in Table [T]!. 
Since we expect the algorithm to converge more quickly in 
this highly idealized case, we decided to continue running the 
GA for 200 generations for the more difficult tests that fol- 
low. The local analysis using LM+SVD brings the final value 
of the metallicity (Z) closer to the source model, while re- 
taining comparable accuracy for the other parameters (see the 
row labeled "local" under H&Hl in Table [T]!. Note that the 
unusually small values of Xr in Table [T] are a consequence of 
adopting scaled-up BiSON errors for the model frequencies, 
which vastly overestimate the true theoretical uncertainties in 
this case. 

3.2. Calibrating Surface Effects 

The biggest challenge to comparing the oscillation frequen- 
cies from theoretical models with those actually observed in 
solar-type stars are the systematic errors due to surface ef- 
fects. The mixing-length parameterization of convection that 
is used in most stellar models is insufficient to describe the 
near-surface layers, and this leads to a systematic difference 
of several ^Hz (up to about 0.3% for a solar model) between 
the observed and calculated frequencies (see Figure |2]l. The 
offset is nearly independent of the spherical degree (I) of the 
mode and grows larger towards the acoustic cutoff frequency. 
The 3D simulations of convection that might in principle re- 
duce this discrepancy for individual stars are far too computa- 
tionally expensive for the model-fitting approach that we are 
developing. Instead, we adopt the me thod for empirical cor - 
rection of surface effects described by lKjeldsen et al.l (l2008h . 
which uses the discrepancies betw een Model S and GOLF 
data for the Sun jLazrek et al.ll 19971) to calibrate the empirical 
surface correction. 

The models that we have adopted for the pipeline include 
a slightly different set of physic al ingredients than what was 
used t o produce Model S from IChristensen-Dalsgaard et al.l 
(11996) "I To characterize the surface effects, we need a model 
that uses the same physics as the other models in the present 
investigation, while matching the Model S frequencies as 
closely as possible. To find such a model (hereafter Model 
S*), we used our pipeline to fit the computed frequencies of 
Model S. Since this involves the comparison of two sets of 
model frequencies — both of which include a mixing-length 
parameterization of convection — there are no surface effects 
to consider. Again we allowed the GA to fit 36 frequencies 
for modes with n = 14-25 and I = 0-2 and we used the same 
scaled BiSON errors and the same constraints on the effective 
temperature and luminosity (see 33.1b . The resulting optimal 
models from the global search and the local analysis had the 
parameter values listed in the rows be low "Model S*" in Ta- 
ble[T] Following lKTeldsen et all dlOOSl) . we fit a power law to 
the differences between the frequencies of the radial modes 
of this model and the corresponding frequencies from BiSON 
data to characterize the surface effects (see Figure |2] and Ap- 
pendix |A|. We found a power law expone nt b = 4.82, sligh tly 
lower than the value (b = 4.90) derived bv lKieldsen et al.l us- 
ing data from the GOLF experiment. With this exponent 
fixed, the recipe of Kjeldsen et al. describes how to predict 
the surface effect for any other set of calculated oscillation 
data, allowing us to apply this empirical correction to each of 

Model S also included diffusion and settling of heavy elements, and was 
based on an eailier version of the OPAL opacities. 
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Fig. 2. — The offset due to surface effects between the Model S* frequen- 
cies and the con'esponding radial modes (circles) from the BiSON data, along 
with the power law fit to these modes (dashed line). The resulting empiri- 
cal correction is also applied to the dipole (triangles) and quadrupole modes 
(squares), since the offset is not a strong function of spherical degree. 



our models before comparing them to observations. 

If our strategy of making this empirical correction to each 
of our models is to succeed, it must not only work well for 
models in a certain region of the search space — it must work 
best for the model that simultaneously matches all of the in- 
dependent observational constraints within their uncertainties. 
Our second H&H exercise (with parameters listed in the row 
labeled "H&H2" in Table ^ was designed to test the behav- 
ior of the pipeline with the surface correction included. The 
source model and other constraints were identical to those 
used for H&Hl, but the surface effect was first calculated 
from Eq. dAlb and then applied to each model frequency prior 
to fitting so that the H&H2 data would mimic real observa- 
tions. We fixed the value of b to 4.82 (as determined above) 
and we calculated the value of a using Eq. dAlb with the scal- 
ing factor r determined using BiSON data as the reference 
set of observations. The fitting procedure then included the 
empirical surface correction — attempting to remove the sys- 
tematic frequency errors that we artificially introduced — as 
detailed in Appendix lAl Again this test reveals how long we 
must run the GA for it to converge to the global solution, but 
it also quantifies any systematic errors on the derived parame- 
ter values that arise from our implementation of the empirical 
surface correction. 

As expected, the GA converges to the global solution more 
slowly with the inclusion of the surface correction. Within 
150-200 generations, two of the four runs converged to the 
parameters listed in the row labeled "global" under H&H2 in 
Table [T] while the other two remained in a nearby local min- 
imum. If we had continued the search, these two runs would 
eventually also have found the global solution. The largest 
systematic error appears in the value of Z, which is about 6% 
(two sample points) higher in the global solution than in the 
source model. There are also smaller errors (^1%) on the 
derived values of Yq and t. The results of the local analy- 
sis (listed in the row labeled "local" under H&H2 in Table [ij 
bring both Z and r closer to the source values while preserv- 
ing the accuracy of the other parameters. The remaining er- 
rors are not a deficiency in the ability of the method to find the 
true solution — the identified model actually results in a lower 
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than the source parameters. 
The reason for this counter-intuitive result is that we can 
only generate a surface correction when there is a reference 
set of observed oscillation frequencies (see AppendixlAli. Our 
source model used the BiSON data for reference, and these 
data have a slightly different value of the large frequency 
separation ({Avq) = 134.673 /iHz) than the resulting source 
model ((Afo) = 134.759 /iHz). Since this quantity is used by 
the binary decision tree to fit the age of the model, the small 
difference leads to a slightly offset derived age — which sub- 
sequently modifies the optimal composition through intrinsic 
parameter correlations. In other words, when generating a 
model with the surface correction included from a given set of 
parameters, those same parameters will generally not provide 
the best fit to the resulting observables. This actually high- 
lights our inability to generate realistic data with surface ef- 
fects, rather than an inherent limitation in our fitting method. 
In any case, this exercise demonstrates that the systematic er- 
rors resulting from our surface correction are small, and that 
running the GA for 200 generations should be sufficient for 
real observations. 

3.3. Validation with Solar Data 

Ultimately, our model-fitting pipeline can only be judged a 
success if it leads to accurate estimates of the stellar properties 
for the star that we know best: the Sun. Up to this point, we 
have essentially been fitting models to synthetic data — using 
the solar data from BiSON only to calibrate the empirical sur- 
face correction and to provide realistic errors. There are many 
other ingredients in our models that could in principle be in- 
sufficient descriptions of the actual conditions inside of real 
stars — deficiencies that could easily lead to systematic errors 
in our determinations of the optimal model parameters for a 
given set of oscillation data. For example, we initially tried to 
use models that employed the simpler EFF equation of state 
(^iggleton et al. 1973) for computational expediency, but this 
led to estimates of the stellar mass about 10% too high for the 
Sun, and unacceptably large systematic errors on many of the 
other stellar properties. Even attempting to ignore the effects 
of helium settling proved to be too coarse an approximation, 
leading to 5% errors on the mass. The only potential ingredi- 
ent that we omitted without serious consequences was heavy 
element settling. This is not to say that simpler stellar models 
cannot be used in the analysis of asteroseismic data, but rather 
that some of the more sophisticated ingredients are required to 
obtain accurate results from a global search of the parameter 
space. 

Having demonstrated the effectiveness of the method by fit- 
ting our models to synthetic data, and after calibrating the 
empirical surface correction using the differences between 
Model S* and the BiSON data, we finally applied our model- 
fitting pipeline to solar data from the BiSON and GOLF ex- 
periments. The oscillation frequencies from these two sources 
are identical to each other within the observational uncertain- 
ties, but their noise properties are slightly different — allowing 
us to quantify any systematic errors that might arise from sub- 
tle effects in the data acquisition and analysis methods. In 
both cases we used the same set of modes referenced in the 
earlier experiments (n = 14-25, I = 0-2) with the respective 
errors again scaled up by a factor of 10, and the same con- 
straints on the effective temperature and luminosity (see 33.1b . 
The two sets of input data differed only in the absolute val- 
ues of the oscillation frequencies (yielding distinct values of 
(Ai'o) for fitting the stellar age), and in the statistical uncer- 
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Fig. 3. — An echelle diagram for the BiSON data (solid points), wliere 
we divide the oscillation spectrum into segments of length (Ai/q) and plot 
them against the oscillation frequency, along with the optimal model from 
our asteroseismic modeling pipeline (open points). Note that the pipeline 
only used the / = 0—2 frequencies between the dashed lines for the fit, but the 
resulting optimal model also matches the 1=3 modes and frequencies outside 
of the fitting range. 



tainties assigned to each mode (leading to subtle differences 
in the weighting of the fit). 

In general, the four independent runs of the GA for each 
data set led to slightly different results after 200 generations, 
so we list only the solution leading to the best model in the 
rows below "BiSON" and "GOLF" in Table[T] We found that 
if we continued running the GA for up to 300 generations, the 
same optimal solution was identified in 2-3 of the indepen- 
dent runs. To reduce the total computing time required for 
future experiments, we decided to stop the GA earlier when 
1-2 of the runs would still reliably identify the global solu- 
tion. Both data sets lead to identical values of the mass and 
metallicity from the global search, with slight variations in 
the values of the other parameters. These minor differences 
largely disappear after the local analysis. Note that because 
we multiplied the true observational errors by a factor of 10 
for the fitting, the resulting values of Xr ^0. 1 . Although 
the fits used a limited range of frequencies and did not include 
I = 3 modes, the optimal models also match the modes with 
lower frequencies and higher degree (see the BiSON fit in Fig- 
ureO and accurately reproduce the known solar properties. 

4. DISCUSSION 

With the successful validation of our model-fitting pipeline 
using solar data, we now need to ensure that our adopted 
treatment of surface effects yields reasonable optimal models 
when applied to other stellar data. While numerous ground- 
based and space-based data sets on other stars are already 
available, the computation-intensive nature of our method 
necessarily limits the present investigation to the single case 
of the Sun. The obvious next step is to use archival ground- 
based data on several well studied solar-type stars to validate 
the pipeline for various stellar masses (e.g. a Cen A & B, near 
^1.1 and ^0.9 M© respectively) and for different evolution- 
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Table 1 . Results of Model-fitting Experiments 
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ary stages (e.g . the "future Sun" /? Hyi, at '^7 Gyr). Although 
iKieldsen et al.l ( [2008,,) demonstrated their method by applying 
it to models of these three stars in addition to the Sun, our 
experience using it with solar data suggests that a global ex- 
ploration of the models may present additional challenges. 

Since our stellar models and the empirical correction for 
surface effects have both been calibrated using a main- 
sequence star at 1.0 Mq, the a Cen system will help vali- 
date the models with interior physical conditions that differ 
slightly from those of the Sun. Its proximity and multiple na- 
ture make it an excellent second test of our pipeline, because 
it has very well deter mined properties inc luding stellar radii 
from interferometry (iKervella et al.l l2003h . There are also 
strong constraints on the component m etallicities and effec- 
tive temperatures (iPorto de Mello et"an f2008). while the ini- 
tial composition and age of the two stars are presumably iden- 
tical. In the next phase of this project, we p lan to use the pub- 
lished oscillation fre quenci es of a Cen A (iBouchy & Carried 
120021; 'Bedding et al."2004'; 'Bazot et al."2007') and a Cen B 
ICarrier & Bourban 2003; Kjeldsen et al. 2005) with the addi- 
tional constraints from interferometry, spectroscopy, and the 
binary nature of the system to further validate o ur pipeline and 
the em pirical correction for surface effects. IKieldsen et al.l 
(l2008h successfully applied their recipe to a set of stellar mod- 
els that broadly resemble the components of the a Cen sys- 
tem, so we have good reason to believe that our implementa- 
tion will also succeed — ^but this remains to be demonstrated. 

The G2 subgiant (3 Hyi has long been studied as a "future 
Sun", with an age near 7 Gyr It has been characterized almost 
as extensively as the a Cen system, incl uding recent inter- 
ferometric measurements of its diameter (iNorth et al.ll2007h 
and dual-site asteroseismic observation s that determined it s 
mean density with an accuracy of 0.6% (Bedding et al.ll2067t) . 
These data included the detection of several 1=1 modes 
that deviate from the asymptotic frequency spacing, suggest- 
ing that they are "mixed modes" behaving like g-modes in 



the core and p-modes in the envelope. This is expected for 
evolved stars like /3 Hyi because as they expand and cool the 
p-mode frequencies decrease, while the g-mode frequencies 
increase as the star becomes more centrally condensed. This 
leads to a range of frequencies where these modes can overlap 
and exchange their character, manifested as so-called avoided 
crossings. This behavior changes very quickly with stellar 
age, and propagates from one mode to the next as a star con- 
tinues to evolve. Consequently, the particular mode affected 
yields a very strong cons traint on the age of the star (see 
IChristensen -Dalsgaard 2004). In subsequent work, we plan 
to use the published oscillation data for /3 Hyi (iBeddi ng et al.l 
l2007h along with the constraints from interferometry and 
spectroscopy to validate our pipeline and the empirical treat- 
ment of surface effects for stars that are significantly more 
evolved than the Sun. This will require an automated method 
to recognize mixed modes in the data set and to incorporate 
them into the optimization of stellar age along each track. 

Once we have validated the model-fitting pipeline with ad- 
ditional stars that sample a range of masses and evolution- 
ary stages, we can begin to consider additional observables 
and parameters that are not constrained by currently avail- 
able data sets. High-quality asteroseismic data are soon ex- 
pected from the Kepler mission, spanning sufficiently long 
perio d s of time that the effects of rotation (iGizon & Solankil 
'2003'; Ballot etal. 2006, 2008) and magnetic activity cy- 
cles (Chaplin et al. 2007; Metcalfe et al. 2007) should be de- 
tectable. The Kepler mission is designed to discover Earth- 
sized habitable planets, and our model-fitting pipeline will be 
able to characterize the planet-hosting stars with asteroseis- 
mology. This is essential to convert precise transit photome- 
try into an absolute radius for the planetary body. In addition, 
accurate rotation rates and ages will provide clues about the 
formation and evolution of the planet-hosting systems. The 
determination of accurate stellar properties for a broad array 
of solar-type stars will give us a new window on stellar struc- 
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ture and evolution, and will provide a broader context for our 
understanding of the Sun and our own solar system. We hope 
to facilitate this process by applying our stellar model-fitting 
pipeline to the data that will soon emerge from the Kepler 
mission. 
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APPENDIX 

A. THE EMPIRICAL SURFACE CORRECTION 

For convenience we summarize the analysis of' Kjeldsen et al.l (1200 8). to correct for the near-surface frequency effects. This is 
based on analyzing the large frequency separation Aj/„ ; (cf. §1) and the individual frequencies. Since the offset from incorrect 
modeling of the near-surface layers is not a strong function of /, we can derive the correction by considering only the radial modes 
(I = 0) and then apply it to all of the modes. 

Suppose we have a set of observed frequencies for radial modes i^^^q^K where n is the radial order Let i'^q'''^^ be the frequencies 
of the model that we seek — the one that best describes the parameters and internal structure of the star, but which still fails to 
model correctly the surface layers. In the solar case, this model is defined by the inversion of the frequencies over a broad range 
of available modes. Here the difference between the observed and best model frequencies is found to be well fitted by a power 
law, which has the convenient property of being free of a frequency scale: 

/ (obs)\ 

-ir--r'=4^ , (Al) 



1^0 



(see also lChristensen-Dalsgaard & Goughlll980h . where is a suitably chosen reference frequency, and a and b are parameters 
to be determined. 

In practice we do not know the best model. Thus, we must determine the correction from a reference model, with frequencies 
vIjqK which is assumed to be close to the best model. In that case we have to a good approximation from homology scaling 

-r^^"^ = -fo^ (A2) 
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where the scaling factor r is related to the mean densities of the best and reference models. 
Substituting Eq. (IA2b into Eq. dAlb and differentiating with respect to n gives 



(obs) 



rAiy' 



,(ref) . 
n.O ■ 



ab 



, Xobs) ' 



b-l 



1^0 



(A3) 



from which we finally obtain 



and 



r = ib 





(ref) 


A-r- 


(obs) 










A-it" 





(A4) 



(A5) 



(for details, seel Kieldsen et al.P 2008'). If we know b then we can calculate r using Eq. (IA4l l. or vice versa using Eq. jASI) . We can 
then obtain a using equations dAll ) and (IA2I) . 

Applying thi s method to the Sun, we us e d the frequencies of Model S,t, which were obtained by fitting the frequencies of 
Model S from Christens en-Dalsgaard et all (119961) . but using slightly modified physics (see 33.21 ). We assumed this to be the 
"best" solar model, whi ch means that w e can set r = 1 (see Eq. IA2b . The observed solar frequencies were obtained from the 
BiSON experiment ( Chaplin et al.lll99^ . From the procedure described above, choosing j/q = 2832.82 /iHz, setting r = 1 and 
using t he data for the do zen radial modes with n = 14-25, we obtained a value of b = 4.82. 

Ki eldsen et alj (l2008h repeated a similar analysis for different numbers of modes (all values from 7 to 1 3), and found the derived 
value of b to range from 4.4 to 5.25. Clearly, the frequency differences do not exactly follow a power law, and so the exponent in 
the power-law fit depends substantially on th e frequency range . However, the value for a varied by less than 0.1 iiHz in all cases. 
If we repeat our analysis using GOLF data (Lazrek et al.'1997) for the same set of modes, we derive a value of b = 4.97. Both of 
our values are comparable to the value = 4.90 derived by Kieldsen et al. (2008). 

After calibrating the value of the exponent b using solar data, the prescription of iKieldsen et al.l describes how to predict the 
surface effect for any other stellar model being compared to real observations. The value of r is determined by applying Eq.( IA4b 
to each trial model through a comparison of the calculated frequencies with the observations. The value of a is then determined 
using equations (lAlb and (IA2I) . and subsequently used in conjunction with b to predict the surface effect for each frequency of 
the model. This recipe is applied to every model calculated during the global search by the GA, and during the local analysis 
using LM+SVD. Effectively we are adjusting the predicted surface effect for each model based on the difference between the 
mean densities implied by the two sets of frequencies. In addition, the magnitude of the surface effect reduces the weight of the 
mode in the final fit, giving more weight to low frequencies but still attempting to optimize the match at high frequencies. The 
final weight for each mode is determined by the quadratic sum of the (statistical) observational error and half the (systematic) 
surface correction. 



B. SINGULAR VALUE DECOMPOSITION 

As discussed in detail by iBrown et al.l (1 19941) . the analysis of the linearized relation between observables and parameters is 
conveniently carried out in terms of Singular Value Decomposition (SVD) because it has some very useful properties when 
applied to a matrix. Suppose we have a model that when given a set of input parameters P (we normally do not know P and 
we wish to retrieve them) produces an output set of expected model observables. Generally the model- (or parameter-) fitting 
problem consists of matching a set of observations O to the output model observables, by making parameter changes ^P to P, 
until we find the model that fits the observations best. Once we find P, to exploit SVD we calculate the partial derivatives of 
each of the observables with respect to each of the model parameters P, and subsequently divide by the measurement errors. The 
result is what we call the design matrix D. 

The SVD of D (= UWV^) neatly describes the relationship between the model observables and the stellar parameters through a 
set of linear transformation vectors. Here the matrices U and V are of order M xN and N xN respectively, where M is the number 
of observables and is the number of parameters; W is a diagonal matrix of order x A^, with diagonal elements consisting of 
the singular values Wk- U and V consist of orthonormal vectors that span the observable and parameter spaces respectively, while 
the diagonal elements of W assign a relative importance to each vector. Each column vector of U describes the content of each 
observable, the singular value Wk quantifies the importance of this vector (i.e., it organizes the observables by amount and type of 
information), while each column vector of V describes how the information from each observable vector is distributed among the 
parameters (i.e., which observables are important for constraining certain parameters). We will see this more clearly below. If the 
models were linear then the solution <5P, in a least-squares sense, could be found simply from the discrepancies SO between the 
observations and the model observables, normalized by the errors. In terms of SVD, the result can be written as ^P = VW"'U^ 
60, which shows clearly that the parameter changes required to reconcile P with the observables comes from a product of (1) the 
discrepancies between the observations and the observables, and (2) a set of linear transformation vectors that are scaled by the 
inverse of the singular values. Below we will explore how we can use SVD to illucidate both the information contained in the set 
of observables and the available constraints on the parameters. 

Figure IB4I illustrates part of the orthonormal matrix U (left) and the full orthonormal matrix V. Because we have N = 5 
parameters in the problem, each of the matrices consists of 5 column vectors; thus we plot and V^, ordered such that the 
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Fig. B4. — The and orthonormal matrices from the SVD of the design matrix, evaluated for the model that best fits the BiSON data. Each row in 
shows how the set of observables are responsible for each of the parameters given in each row of V^, e.g., row 1 contains mostly the parameter of Z, and this 
is constrained mainly by the 3 frequencies that appear in row 1 of U^. If we scale the rows in V'^ by the inverse of the singular values we obtain the parameter 
correlations, which describe how to define the 5-dimensional en'or boundary. 



horizontal vectors corresponding to the largest singular values Wk are at the top. Each column vector of U has M = 38 elements, 
which is the number of observations that are used to determine the solution, but we only include the first 5 elements for illustration. 
The elements of the matrix can only take values between -1 and +1, and each number is represented by a triangle whose size and 
direction are proportional to the magnitude and sign of each element. The column vectors are labeled 1-5: 1 represents the first 
column vector, which is associated with the largest singular value, 2 is the second largest singular value, while 5 is the smallest. 
Thus, the vectors are effectively ranked according to their information content, where the least important information is given in 
the vectors corresponding to the smallest Wk- For example, if we were to perform the matrix multiplication of UWV^ using only 
the top vectors, we would be able to recover most of the input matrix D. Essentially this is equivalent to zeroing the values of 
W that are below a certain threshold before doing the multiplication. Clearly the observables that appear among the top vectors 
are most important for constraining the parameter solution, while the parameters that appear among the top vectors are the most 
constrained. 

Inspecting the left panel of Figure IB4I by scanning down the fourth column we see that the observable Teff has components 
mainly in the second and fourth vectors. The luminosity is present in vectors 2-4, while the individual frequencies are present 
mostly in vectors 1 and 2. The most important information comes from the first vector, and we can see clearly that in this case the 
seismic information dominates. We also see that and Tetf usually appear in the same vectors, implying that these observables 
work together to constrain the parameter solution. The right panel shows how the information from the observables is distributed 
among the parameters. The top vector has mainly one large component, Z. This implies that Z is most constrained by the set of 
observables, and in particular by the observables that appear in the top row of the matrix, i.e. the frequencies. Similarly, the 
mixing-length parameter (a) appears predominantly in the fourth vector, and is constrained mainly by Tgff and L*. 

Returning to the linear parameter-fitting problem, we saw that calculating depended strongly on the size of the singular 
value associated with the transformation vectors (as well as the discrepancies). Because Z appears in the first vector and the 
inverse of the first singular value is the smallest, this means that Z is not allowed to deviate much from its initial value. So 
if we find that we are in a local minimum that is not the global solution, it is unlikely that we will find a better solution that 
deviates much in Z. For this reason we choose to explore a range of Z values during the local analysis. The vectors Vr are the 
basis vectors, and describe the normalized changes that we can make to each of the parameters while respecting the constraints 
imposed by the observables. The inverse of the Wk describe the magnitude of each of the vectors that keeps the parameters within 
1(7 of the final answer Since the vectors corresponding to larger k have smaller Wk (larger inverses), the parameters associated 
with these vectors will extend further into the parameter space, i.e. they are less constrained by the observations. 

Inspecting the matrix in Figure IB4I note that vector 2 reveals both and Yo to be positive, indicating that they should 
be increased or decreased together to satisfy the constraints imposed by the observables. However, vector 3 implies that if 
increases, Yo should decrease. When both of these vectors are projected onto the M^-Kq plane and scaled according to I/W2 and 
I/W3, vector 3 is the longer of the two. These scaled vectors are shown as arrows in Figure lB5l defining the semi-major and 
minor axes of the projected la error ellipse. The error ellipse defines the boundary of the allowed parameter values that will 
change the value by less than Ax^, where Ax^ is the value that defines the la error for the relevant number of degrees of 
freedom. This figure shows clearly that there is a very strong inverse co rrela tion between these parameters, and hints at the size 
of the uncertainty in each one". Though not as clearly visible in Figure lB4l Z and Fq have a strong direct correlation, so for the 
local analysis we explore a range of Yo as well as Z. 



' ' The correlation arises naturally from the fact that increasing either M^, or Yq tends to increase the luminosity, so if one of these parameters is increased the 
other must be decreased to satisfy the luminosity contraint while continuing to match the other observables. 
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Fig. B5. — The longest projections of the 5-dimensional error boundary on the 2-dimensional parameter space of and Yq, as described by SVD, for the final 
fit model to the BiSO N da ta. The arrows are the 2-element directional vectors given by the projections of V2/W2 and V3/W3 which dominate the determination 
of Mt and yb (cf. Fig. lB4l . and the ellipse is defined by these vectors. 



Finally, the covariance matrix can be expressed in a very compact form using SVD; the covariance between parameters Pj and 
Pi is: 



In particular, the formal uncertainties criPj) on each of the parameters are obtained from <j(Pj)^ = J2k^fk/^kk- Given this very 
clear form for describing the correlations between the parameters and their uncertainties, we can also determine the uncertainties 
on the model observables (such as and L*) or some other properties (such as the core hydrogen content) by calculating the 
models at the boundaries of the error ellipse. Equation ( IB lb can also be used to investigate how to reduce the uncertainties in 
each of the parameters — with what precision do we need to measure the individual observables? What other information may be 
useful to help break the degeneracies between some parameters? And are there limits to what we can learn from the available 
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